function [LE1, LE2] = LEs(func, x0, q0, a, b, c, k, N, transient)
% LEs 计算二维映射的两个 Lyapunov 指数，模型迭代函数由 func 指定
%
% 输入：
%   func      - 模型迭代函数句柄，其调用格式为 [x_next, q_next] = func(x, q, a, b, c, k)
%   x0, q0    - 初始状态
%   a, b, c, k- 模型参数
%   N         - 用于计算 LE 的有效迭代步数（舍去暂态后的步数）
%   transient - 舍去的暂态步数
%
% 输出：
%   LE1, LE2  - 计算得到的两个 Lyapunov 指数

    % 初始化状态与扰动矩阵 Q（2x2 单位阵）
    x = x0;
    q = q0;
    Q = eye(2);
    sum_log = zeros(2,1);
    
    total_steps = transient + N;
    for i = 1:total_steps
        % 迭代系统：调用传入的 func 进行更新
        [x, q] = func(x, q, a, b, c, k);
        % 计算当前雅可比矩阵
        J = jacobian(x, q, a, b, c, k);
        % 将扰动向量传播
        Q = J * Q;
        % QR 分解，对 Q 进行正交化
        [Q, R] = qr(Q);
        
        % 舍去暂态后累计对角线元素对数
        if i > transient
            for j = 1:2
                sum_log(j) = sum_log(j) + log(abs(R(j,j)));
            end
        end
    end
    
    % 取平均得到 Lyapunov 指数
    LE = sum_log / N;
    LE1 = LE(1);
    LE2 = LE(2);
end
